In order to plot fields, we run the following commands:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
We import the math library:
In [2]:
import math
The modules to import for dealing with grids are:
In [3]:
from pygsf.mathematics.arrays import *
from pygsf.spatial.rasters.geotransform import *
from pygsf.spatial.rasters.fields import *
For calculating pathlines:
In [4]:
from pygsf.space_time.movements import interpolate_rkf
The example circular motion vector field has components:
v = y i - x j
as deriving from the equation:
v = - z x r
where z is the vertical vector, r the position vector and x the vector product.
In [5]:
k = 2 * math.pi
def z_func_fx(x, y):
return k*y
def z_func_fy(x, y):
return -k*x
The velocity field parameters for testing the results are:
v = w * r
w = v / r
|v| = sqrt(k^2y^2 + k^2x^2) = k * r
1 cycle -> 2 pi r
v = ds / dt -> ds = v * dt
2 pi r = v dt
2 pi r = v T -> T = 2 pi r / v = 2 pi / k
In [6]:
rows=100; cols=100
In [7]:
size_x = 1; size_y = 1
In [8]:
tlx = -50.0; tly = 50.0
In [9]:
gt1 = GeoTransform(
inTopLeftX=tlx,
inTopLeftY=tly,
inPixWidth=size_x,
inPixHeight=size_y)
In [10]:
fx1 = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_fx)
In [11]:
print(fx1)
In [12]:
fy1 = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_fy)
In [13]:
print(fy1)
To visualize the parameters of the flow, we calculate the geographic coordinates:
In [14]:
X, Y = gtToxyCellCenters(
gt=gt1,
num_rows=rows,
num_cols=cols)
and the vector field magnitude:
In [15]:
magn = magnitude(
fld_x=fx1,
fld_y=fy1)
In [16]:
fig = plt.figure(figsize=(14, 6))
plt.contourf(X, Y, np.log10(magn), cmap="bwr")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Magnitude (log10)')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.axis("image")
plt.title('Vector field magnitude and streamlines')
plt.xlabel('x')
plt.ylabel('y')
Out[16]:
In [17]:
import math
from pygsf.spatial.vectorial.vectorial import Point
from pygsf.spatial.rasters.geoarray import GeoArray
In [18]:
ga = GeoArray(
inGeotransform=gt1,
inProjection="undef",
inLevels=[fx1, fy1])
In [19]:
time_increm = 1.0e-4
In [20]:
period = 2 * math.pi / k
In [21]:
number_of_cycles = 100
In [22]:
steps = number_of_cycles * (period / time_increm)
print (steps)
In [23]:
first_pt = Point(0, 20)
In [24]:
str_pt = first_pt
pts_x, pts_y = [first_pt.x], [first_pt.y]
for n in range(int(steps)):
end_pt, error = interpolate_rkf(
geoarray=ga,
delta_time=time_increm,
start_pt=str_pt)
if end_pt is None:
break
pts_x.append(end_pt.x)
pts_y.append(end_pt.y)
str_pt = end_pt
print (end_pt)
After 100 cycles the calculated point position is in the expected (initial) position: x=0, y=200.
In [25]:
fig = plt.figure(figsize=(14, 6))
plt.contourf(X, Y, np.log10(magn), cmap="bwr")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Magnitude (log10)')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.scatter(pts_x, pts_y, color='b')
plt.axis("image")
plt.title('Vector field magnitude and streamlines')
plt.xlabel('x')
plt.ylabel('y')
Out[25]: